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i ABSTRACT 

> ' 

on : 

In this work, we first establish a simple procedure to obtain with 11-figure 

00 ■ 

accuracy the values of Chandrasekhar's H-f unction for isotropic scattering using 
<35 a closed-form integral representation and the Gauss- Legendre quadrature. Based 

on the numerical values of the function produced by this method for various 
7— i . combinations of vuo, the single scattering albedo, and fi, the cosine of the azimuth 

J> . angle 9 of the direction of radiation emergent from or incident upon a semi-infinite 

•i-H . 

^ , scattering-absorbing medium, we propose a rational approximation formula with 

/i 1 / 4 and y/1 — xxj as the independent variables. This allows us to reproduce the 
correct values of H{wq,h) within a relative error of 2.1 x 1CT 5 without recourse 
to any iterative procedure or root-finding process. 

Subject headings: radiative transfer — iJ-function — approximations 



1. Introduction 



Chandrasekhar's H-functions H(wQ,fi) are used to express the emergent intensities of 
radiation reflected by semi-infinite, homogeneous media. They satisfy the following nonlinear 
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integral equation (see, e.g., IChandrasekhar Ill960l ): 



f 1 ^(v) 

H(w , fjt) = l + fJ,H(zu , fj) / — ■ — H(w , r])drj 

Jo V + V 



X 



where zu is the single scattering albedo, and which is also a function of zu , is the char- 
acteristic function specific to the type of scattering the radiation undergoes. We developed, 
for instance, a compact computational method to find the numerical solutions of Eq. (pQ) for 
some representative types of anisotropic scattering in terms of the roots of the characteristic 
equation involving values of \l/ evaluated at the q uadrature points and the associated weights 
of integration employed (IKawabata et al. 1119911 ) . 



In practical applications of the theory of radiative transfer, however, we often need a 
quick and yet relatively accurate approximation for the if-function, especially for isotropic 
scattering. In the present study, we therefore focus on developing a simple numerical proce- 
dure to calculate accurately (at least to the 11th digit) the values of H{wq, //) for conservative 
as well as non-conservative isotropic scattering. We will then construct a rational approxima- 
tion formula based on fitting to the reference values of H(wq, fi) produced with the foregoing 
method. This permits us to avoid the iterative methods that are often employed in this type 
of computations. 

We must note that several usef ul approximation formula e have been proposed by var- 
ious authors. The formula given by iKaranjai and Sen I (119711 ) for arbitrary values of vjq is 
sufficiently compact, but its use requires the solution of a transcendental equation fo r a given 



value of vjq, which we want to av oid. The same is true for the formula derived by iDomke 
(119881 ). The formula proposed by IHapke I (119931) is also handy, but its accuracy is around 
7.7 x 10~ 3 . The formula developed by lKaranjai and Karanjai I (Il99ll ) is sufficiently accurate, 



with a maximum relative error of about 10 -4 , and requi res no root-finding proc ess, but is 
valid only for conservative scattering. Recently, however, iDavidovic et al. I (120081 ) obtained 
a new formula for an arbitrary set of (wq,/j,), which is an order of magnitude more accu- 
rate (the maximum relative error is 7 x 10 -4 ) than any of these. Therefore, it must be an 
interesting challenge to find a formula that is at least an order of magnitude better than 
this. We will make all the numerical computations required for this objective exclusively in 
double-precision arithmetic. 
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2. Formulations and Numerical Computations 
2.1. Integral Representation 



H(w , ji) = exp 



The H-function for iso tropic scattering H(zuq, h ) can be expressed in a closed-form 
integral representation (see iRutily and Bergeat 1119871 ): 

- / In (1 - O7 £ arccot £) ■ ^ 
71 Jo V + 6 

Substituting arccot^ = x in Eq.(|2J), we have 

r'\ (i + cot 2 x) 

/ m(l — wqx cot xj • -— 5 — 

/ (/r + cor x 



H{vjq, ji) = exp 



-dx 



(2) 



(3) 



Eq.flS]) is advantageous for numerical computations because the domain of integration 
is finite. We would therefore like to use this expression along with the Gauss-Legendre 
quadrature to generate the reference values of H(tuq,h) for the purpose of comparison. For 
non- conservative scattering, for which wq < 1, the numerical integration inside Eq.([3]) can 
be carried out without difficulty by means of the quadrature, and the resulting values of 
H(wo, /jl) should be accurate enough even if a single quadrature is applied to the entire do- 
main of integration [0, vr/2]. However, conservative scattering with vjq = 1 poses a numerical 
problem in that the factor 1 — x cot x, inside the natural logarithm involved in the integrand, 
diverges as x tends to 0, which could significantly degrade the numerical accuracy of the 
resulting H-function. With this in mind, for conservative scattering, we divide the domain 
of integration into two parts, [0, e] and [e, 7r/2], where £< 1: 



tt/2 



ln(l — x cot x) 



(1 + cot 2 x) 
(fi 2 + cot 2 x) 



dx — Ji + J 2 



where 



and 



h = ln(l — x cot x) 
Jo 



(1 + cot 2 x) 
(/i 2 + cot 2 X 



■dx, 



(4) 



(5) 



tt/2 



I s 

N 

E 



ln(l — x cot x) 
ln(l — Xj cot Xj) 



(1 + cot 2 x) 
(/i 2 + cot 2 x 

(1 + COt 2 Xj) 
(fl 2 + COt 2 Xj 



■dx 



(e < Xj < -), 



(6) 



where Xj and Wj are the j-th quadrature point for evaluating the integrand and the cor- 
responding weight of the Gauss-Legendre quadrature applied to the x- interval of [e,n/2], 
respectively. 
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Since < x < e < 1 for Ji, we may expand its integrand in a series in x. Retaining the 
terms up to and including the order of x 3 , we get 



ln(l — x cot x) 



1 + cot 2 x 
fi 2 + cot 2 x 



In h 

3 



I5 + * 1 



2 Mn- 



x 2 + o[xf 



(7) 



Substitution of the right-hand side of Eq.flT]) into Eq.(JS]) yields 



. ar 

In h 

3 



- + (l-/i 2 )ln- 
15 K * ' 3 



x 



dx 



£ 

45 
+o[ef 



{30(3 + Ae 2 )ln £ + [1 - 5A • (2 + 31n3)]e 2 - 45(2 + ln3)} 



(8) 



with A 



(see also Appendix A for a higher-order approximation). Because we are 
interested in producing the numerical values of H(l, fi) with 11-figure accuracy, and we have 
the constraint e < 10 -2 , let us set e = 10 -3 for wq = 1 and e = otherwise. The value of 
H(wo,fi) owing to the integral representation is then given by 



H- mteg (w , fx) = exp 



7T 



(9) 



Note that a similar integral representation was employed by iDavidovic et al. I ( 120081 ) to 
produce the reference numerical values of H{wq ) li). However, their domain of integration 
was [0, oo], which makes it necessary to introduce a special scheme of integration. In view 
of this, we believe that the present technique is much more straightforward and easier to 
handle 0. 



2.2. Computational Results with Integral Representation 



We varied the degree iV of the Gauss-Legendre quadrature to calculate the values of 
H{w n , /x) for all the co mbinations of 48 values of zu (0.01 plus 47 values employed by 
( 2008h for their Table 1) and 21 values of \i (0 plus 20 values employed by 
(120081 ) for their Table 1). Note that similar numerical tables are also given 



Davidovic et al. 



Davidovic et al. 



1 The an onymous referee kind ly directed our attention to an alternative integral representation (see, e.g., 
Eq.(16) of iDas and Bera 1 120071 ) . Based on some numerical tests of this formula, we were however led to 
the conclusion that the use of Eq.(3) coupled with Eq.(4) is approximately two orders of magnitude faster to 
produce the values of H(l,fj.) that are correct to the 10th decimal place. A brief discussion on this matter 
is given in Appendix B. 
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by iBosma and de Rooij I (119831 ). but for the combinations of roo 
and 1 and for /x from to 1 with a step of 0.1. 



0.5, 0.7, 0.9, 0.99, 0.999, 



We found that N = 100 is sufficient to br ing o ur results for H(m n , /x ) into complete 
agreement with those of iDavidovic et al. I (120081 ) and IBosma and de Rooij I ( 119831 ) down to 
the 10th decimal place (11 figures altogether). On the other hand, even if the 300th-degree 
Gauss-Legendre quadrature is employed for the single interval [0, 7r/2], the accuracy of the 
resulting value of the H-function for the conservative scattering is much lower: we obtain, 
e.g., H(l, 1) = 2.9077901976 instead of the reference value 2.9078105291. This indicates the 
effectiveness of Eq.Q. 

Next, to further assess the quality of the present approximation, we examined do, the 
zeroth moment of the if-function, which can be expressed in terms of the single scattering 
albedo zu as 

2 



r 1 2 

a = H(xu ,fjb)dfj,= — [1 - vT- ro ol 
Jo ^0 



(10) 



for isotropic scattering (IChandrasekhar Ill960l ). Using the values of H(xuq,h) given by our 
present method, we performed the numerical integrations required to obtain a , again using 
the A^GL-th degree Gauss-Legendre quadrature: 



a 



0; fJ-j) 



w 



Unfortunately, however, with A" = 100, which was sufficient for computing H(vj$, /x), we 
found it impossible to get the value of ao correct to the 9th decimal place regardless of how we 
chose the value of Aql- For instance, the best value we could obtain was a = 2.0000000019 
with Aql = 130 for wq — 1. To obtain values of «o correct down to the 9th decimal 
place, thereby allowing for at most a unit difference at the 10th decimal place relative to 
the theoretical values 2(1 — y/1 — vjq) /zjq, we found it necessary to employ A^ = 300 and 
Aql > 270. Allowing for an adequate margin, we therefore adopted A^gl = 350, with which 
we finally obtained most of the computed values of «o correct to the 10th decimal place, 
e.g., a = 2.0000000000 for w = 1 and 1.4854314511 for w = 0.88. The only exception 
was that for w = 0.9996, 0.9995, 0.999, 0.995, 0.993, 0.965, 0.95, 0.93, and 0.75, the figures 
at the 10th decimal place were larger than the theoretical values by unity. We nevertheless 
concluded that our primary objective was accomplished to a sufficient degree at this stage. 
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2.3. Approximation Formulae 

Now that we have established a highly reliable means of obtaining the numerical val- 
ues of H (mo, fi) , we proceed to the next stage and seek a fast and yet reasonably accurate 
approximation formula for H(zu ,fi). For this, we shall proceed in two steps: (1) try to 
construct a polynomial approximation formula for H(l,fi) accurate to at least five figures, 
and (2) develop a rational approximation formula for H(wq,h) by using that obtained in 
(I)- 



2.3.1. Approximate Formula for Conservative Scattering 

Using polynomials of various degrees K x of /i 1//4 , we made least square fittings to the 
reference values of f7i nteg (l, \x) tabulated for 501 equally spaced values of /i between and 1: 

Ki 

H iQteg (l,fi) = J2 A ^ k , (12) 

k=0 

where x = /i 1//4 , as mentioned above, and are the constants to be determined. The 
choice of the independen t variable x stems from the experience gained in our foregoing work 



(IKawabata et al. Ill99ll ). 



P 



After some experimentation, we found that K\ — 8 yields a satisfactory fit to the 
reference values. The polynomial approximation formula thus obtained is as follows. 

# ap p(l,/f) = 9.999982706853756 x lO^ 1 + 3.465443224211651 x 10" 4 x 

-1.411107006687451 x 10~V + 3.269177042230116 x lO^x 3 
4.133809356648527x 4 - 7.188546622876579x 5 + 7. 772939980710241a; 6 
-3.883055730606847x 7 + 7.595128286312914 x lO" 1 ^ 8 . (13) 

A comparison between the values of i7 app (l,//) generated by t he present formula and 



those obtained with the integral representation H- mieg (l,fi) (see also iDavidovic et al. 112008 
Table 1) indicates that the maximum relative error is 2 x 10~ 6 . The approximate values 
are found to be correct to at least the fifth decimal place, and even the figure at the sixth 
decimal place differs from the correct one by no more than one unit. 
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2.3.2. Approximate Formula for N on- Conservative Scattering 

For a set of L equally spaced values of tuq, we compute the values of if app (l, /x) / Hi Qteg (zuo, /x) 
for M values of /x, which are taken to coincide with the Chebyshe v collocation p oints to apply 
the Chebyshev polynomial approximation method described in iPress et al. ] Jl992h . Then, 



they are approximated by a K 2 -th degree polynomial of x(= /x 1//4 ), as before: 

# apP (l, fj) I H integ (zu , /x) = ) j C k (-cj )x k . (14) 

k=0 

For a given value of wo, the values of the coefficients C k (wo) (k = 0, • • • ,K 2 ) are 
determined by fitting the right-hand side of Eq. (fH]) to the values of if app (l, /x) / H integ (zu , /x) 
computed at M discrete points of /x, as indicated above. As a result, we obtain K 2 sets of L 
values of C k {wo). 

For actual computations, we adopted M = 3, 500 and L = 10, 001 so that the value of 
too ran from zero to unity with a step size of 10~ 4 . For each k, the 10,001 values of Ck{zuo) 
obtained with the Chebyshev polynomial approximations were finally approximated by an 
iVp-th degree polynomial of y/1 — w : 

C k (w ) = J2B k ,nV n , (15) 

n=0 

where 77 = a/1 — w and the coefficients B k n were determined by the standard least square 
method while varying the value of iVp from 4 through 9. The use of y/1 — cc as the inde- 
pendent variable is based on a premise draw n from another past work on the if-function for 



isotropic scattering (jKawabata et al. 1119921 ) 



To assess the quality of the resulting formula, we computed 

K 2 

H app (w , fjt) = # app (l, n)/ ^2 C k (zu )x k (16) 

at the same 48 x 21 grid points on the (wq, /x)-plane as used in Section 2.2, which were then 
compared with those of H integ (wo, /x). The best formula was found with K 2 = 8 and iVp = 8, 
whose coefficients C k (?&o) (k — 1, • • • , K 2 ) are shown below. 



Co = -1.368687418901498 x 10" 6 + 6.744526217097578 x 10~ 5 r/ 

-8.816747094601710 x 10~V + 4.731152489223286 x 10~ 3 r/ 
-1.352739541743824 x 10" 2 r/ 4 + 2.236433018731980 x 10" 2 r/ 
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-2.147081702708310 x 10 _ V + 1.112257595951489 x 10 _ V 
-2.406003988429531 x 10~V (17a) 
Ci = 8.737822937355147 x 10~ 5 - 5.250514244222347 x 10~ 3 r] 
+7.644952859355422 x 10 _ V - 4.664908220536214 x 10~V 
+1.48268819832583977 4 - 2.663033364728811^ 
+2.727252555244034?? 6 - 1.485444888951274r? 7 

+3.340921510758153 x 10~V (17b) 

C 2 = -1.427222952750036 x 10~ 3 + 9.300028322140796 x lO" 2 ^ 
-1.413069914567426r? 2 + 8.880428860986575^ 
-2.866825946137678 x 10+V + 5.178036196746675 x 10+V 
-5.307180734532348 x 10+V + 2.885782084328829 x 10+V 
-6.471219440031649r? 8 (17c) 

C 3 = 9.066801756884433 x 10~ 3 - 6.354984995808299 x lO^r] 

+1.021262226727643 x 10+V - 6.444360574298017 x 10+V 
+2.105330190640368 x 10+V - 3.824039368443171 x 10+V 
+3.930240665640704 x 10+V - 2.139686267143788 x 10+V 
+4.800025272319539 x 10+V (17d) 

C 4 = -2.855922558150419 x 10" 2 + 3.880224653851042^ 

-3.174231079700075 x 10+V + 2.303877926374539 x 10+V 
-7.626655021168267 x 10+V + 1.394034249890738 x 10+V 
-1.438476211044276 x 10+V + 7.852393856327993 x 10+V 
-1.764969590005163 x 10+V (17e) 

C 5 = 4.941209676842531 x 10~ 2 - 3.976393849244121^ 

+6.000178277203062 x 10+V - 4.542543148444882 x 10+V 
+1.512146625692455 x 10+V - 2.779737284749243 x 10+V 
+2.880598698878311 x 10+V - 1.577451021926768 x 10+V 
+3.554375808436865 x 10+V (17f) 

C 6 = -4.798519468590785 x 10~ 2 + 4.112841572654386?? 

-6.655808348671680 x 10+V + 5.000349699512032 x 10+V 
-1.672172432180451 x 10+V + 3.091851778649070 x 10+V 
-3.218110914157008 x 10+V + 1.768094273655673 x 10+V 
-3.994358424590589 x 10+V (17g) 
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C 7 = 2.461700902387896 x 10~ 2 - 2.233648393380449^ 



+3.900465646584139 x 10+V - 2.880699974056035 x 10+V 
+9.688954523412610 x 10+V - 1.802235503900686 x 10+V 
+1.883990440310628 x 10+V - 1.038462482861755 x 10+V 
+2.352061082130820 x 10 + V 



(17h) 



C 8 = -5.211353622987505 x 10" 3 + 4.967427514273564 x lO^r] 



-9.292147966163522?? 2 + 6.773895398390997 x 10+V 



-2.294206635762768 x 10+V + 4.292903843888321 x 10+V 
-4.506396634901928 x 10+V + 2.491623632369491 x 10+V 



-5.657192709351447 x 10+V 



(17i) 



The maximum relative error of the approximate values if app (G7o, /x) given by this formula 
on the grid points of (zu , fi) is 2.1 x 10~ 6 and occurs at tu = 0.996 and \i = 0, whereas 
the mean of the relative errors is 3 x 10~ 7 . The differences between the formula values 
and those of H integ (zu , ji) are within ±4 at the sixth decimal place. If rounded at the fifth 
decimal place, they agree with each other to the last digits except for the three cases with 
(u7 , n) = (0.9, 0.3), (0.995, 0.95), (0.999, 0.7), and (0.9998, 0.4), where the figure in the fourth 
decimal place differs from the correct number by one unit. The maximum relative error of 
the computed values for the zeroth moment is 1.64 x 10~ 7 , arising at wq = 1, where we 
have 2.0000003285 instead of 2 by means of the 350th-degree Gauss-Legendre quadrature. 



We have developed a numerical procedure to evaluate with 11-figure accuracy the value 
of Chandrasekhar's if-function for isotropic scattering for arbitrary sets of the single scat- 
tering albedo, wq, and the cosine of the zenith angle of the emergent or incident direction of 
radiation, /x, using the integral form expression. This should prove useful not only for check- 
ing the accuracy of numerical values of H(zu , /i) computed by other types of approximations 
but also as a practical computational tool in applications. 

The rational approximation formula constructed on the basis of the reference data com- 
puted with the integral-form representation for H(mo,[i) is significantly more accurate than 
any available in the literature, which may well make up for the fact that it is longer than 
others. 



3. 



Conclusion 
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A. Ii up to the e 5 Terms 



Expanding the left-hand side of Eq.flT]) in a series in x and keeping the terms up to and 
including the order of x A , and analytically carrying out the integration of Eq.©, we get 



(2 In £ - 
617 



15750 



2 - In 3)e + — [30(1 - fi 2 ) In e + 5/i 2 (2 + 3 In 3) - 3(3 + 5 In 3)] e z 
+ ^ - ^/x 2 (9 + 25 In 3) + ^//(2 + 5 In 3) 



15 



(2 - 5/i 2 + 3/i 4 ) In e 



e 5 + o[e) 



(Al) 



If £ 



10 6 as in our numerical computations, Ii = —1.691412801800870 x 10 
0, -1.691412671960849 x 10~ 2 for // = 1/2, and -1.691412282441017 x 10" 2 
1 in contrast to -1.691412801800667 x 10" 2 , -1.691412671960754 x 10~ 2 , and 
-1.691412282441016 x 10" 2 respectively obtained from Eq.®. This fact indicates that the 
latter expression for l\ is of sufficient accuracy for our present purposes. 



for \i 
for \i 



B. Alternative Integral Representation for Conservative Scattering 



In the case of zu = 1, we have an interesting alternative to Eq.(3) (IDas and Bera 1120071 ): 



where 



and 



H(l,n) = \/3(l + /i)exp 



9{x) = — atan2 



7T 



1 dx 

e(x)-=- 

X + fJ, 
1 + X 



n 1 1 1 

— x, 1 x In 

2 ' 2 l-x 



atan2 (y, x) 



arctan (y/x) x > 

7r + arctan (y/x) y > 0, x < 

—ii + arctan (y/x) y < 0, x < 

tt/2 y > 0, x = 

-7I-/2 y < 0, x = 

undefined ?/ = 0, x = 



fBll 



(B2) 



(B3) 
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Some remarks on Eq.(Bl) must therefore be in order in comparison with Eq.(3): first 
of all, the integrand 8(x)/(x + fi) involved in Eq.(Bl) exhibits a very steep increase as we 
come close to x — 1 and converges to a limiting value of 1/(1 + fi) at x — 1. Furthermore, 
in the vicinity of x — 0, it also rises almost vertically from to nearly 0.5 over a short range 
in x for small but non-zero values of [i. This fact makes the numerical integration of this 
integrand extremely difficult particularly if /i is small. 



For numerical tests of the efficiency of Eq.f lBip . we have widely varied the degree N of the 
Gauss- Legendre quadrature employed for carrying out the required integration. The resulting 
values of H(fj,) are shown in Table 1 for N = 300, 3000, and 30000 in comparison with 
those obtained from Eq.(3) with N = 300 coupled with Eq.(4) . Also shown in the column 



design ated by " Das-Bera" are the computational values taken from Table-8 of iDas and Bera 



( 120071 ). who used the Simpson's one third rule (the number of division points is unknown). 



It is now obvious that we need to employ the Gauss-Legendre quadrature with iV > 
30000 to get the values of H(n) by means of Eq. (1Blj) correct down to the 10th decimal 
place, and that Eq. (3) together with the procedure discussed in the text of the present work, 
requiring only N = 300 or less , is significantly faster at least for the purpose of generating 
reference numerical values. 



Table 1: Accuracy comparison of the values of H(/j,) for conservative scattering 



/l 


Eq.(3) 


Das-Bera 


N = 300 




N = 3000 


N = 30000 


0.00 


1.0000000000 


1.0000000000 


1.0000000446 


1 


0000000002 


1 


0000000000 


0.05 


1.1365748468 


1.1365748417 


1.1365748951 


1 


1365748471 


1 


1365748468 


0.10 


1.2473504425 


1.2473504371 


1.2473504930 


1 


2473504428 


1 


2473504425 


0.20 


1.4503514128 


1.4503514071 


1.4503514667 


1 


4503514131 


1 


4503514128 


0.30 


1.6425222645 


1.6425222585 


1.6425223208 


1 


6425222648 


1 


6425222645 


0.40 


1.8292756032 


1.8292755970 


1.8292756614 


1 


8292756035 


1 


8292756032 


0.50 


2.0127787700 


2.0127787636 


2.0127788298 


2 


0127787703 


2 


0127787700 


0.60 


2.1941330193 


2.1941330128 


2.1941330805 


2 


1941330197 


2 


1941330193 


0.70 


2.3739749125 


2.3739749059 


2.3739749748 


2 


3739749129 


2 


3739749125 


0.80 


2.5527043168 


2.5527043101 


2.5527043801 


2 


5527043172 


2 


5527043168 


0.90 


2.7305876649 


2.7305876581 


2.7305877289 


2 


7305876652 


2 


7305876649 


1.00 


2.9078105291 


2.9078105222 


2.9078105939 


2 


9078105294 


2 


9078105291 
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